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A  DOMAIN  WALL  MODEL  FOR  HYSTERESIS  IN  PIEZOELECTRIC  MATERIALS 

RALPH  C.  SMITH*  AND  ZOUBEIDA  OUNAIESt 

Abstract.  This  paper  addresses  the  modeling  of  hysteresis  and  nonlinear  constitutive  relations  in 
piezoelectric  materials  at  moderate  to  high  drive  levels.  Hysteresis  and  nonlinearities  are  due  to  the  domain 
structure  inherent  to  the  materials  and  both  aspects  must  be  addressed  to  attain  the  full  potential  of  the 
materials  as  sensors  and  actuators  in  high  performance  applications.  The  model  employed  here  is  based 
on  previously  developed  theory  for  hysteresis  in  general  ferroelectric  materials.  This  theory  is  based  on  the 
quantification  of  the  reversible  and  irreversible  motion  of  domain  walls  pinned  at  inclusions  in  the  material. 
This  yields  an  ODE  model  having  five  parameters.  The  relationship  of  the  parameters  to  physical  attributes 
of  the  materials  is  detailed  and  algorithms  for  determining  estimates  of  the  parameters  using  measured  values 
of  the  coercive  field,  differential  susceptibility  and  saturation  properties  of  the  materials  are  detailed.  The 
accuracy  of  the  model  and  its  capability  for  the  prediction  of  measured  polarization  at  various  drive  levels 
is  illustrated  through  a  comparison  with  experimental  data  from  PZT5A,  PZT5H  and  PZT4  compounds. 
Finally,  the  ODE  model  formulation  is  amenable  to  inversion  which  facilitates  the  construction  of  an  inverse 
compensator  for  linear  control  design. 

Key  words,  hysteresis  model,  piezoelectric  materials 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Piezoelectric  materials  provide  the  capability  for  designing  actuators  and  sensors 
which  are  compact,  lightweight,  can  be  molded  or  constructed  in  a  variety  of  configurations,  and  are  relatively 
inexpensive.  Hence  the  materials  are  being  employed  in  an  increasing  number  of  structural  and  structural 
acoustic  applications  with  uses  which  include  active  vibration  control,  the  attenuation  of  structure-borne 
noise,  micropositioning,  and  high  performance  structural  drivers.  The  mechanisms  which  provide  the  ma¬ 
terials  with  both  their  sensor  and  actuator  capabilities  are  due  to  the  noncentrosymmetric  nature  of  the 
materials  and,  more  specifically,  to  domain  switching  in  response  to  applied  fields  or  stresses.  In  the  former 
case,  the  polar  changes,  which  occur  when  ions  displace  to  align  with  an  applied  field,  produce  the  strains 
used  to  actuate  the  underlying  structure.  Conversely,  the  application  of  stresses  produces  deformations  in  the 
material  which  alter  the  polarization  and  subsequently  generate  the  voltages  measured  when  the  materials 
are  employed  as  sensors.  These  are  the  converse  and  direct  piezoelectric  effects. 

As  a  result  of  the  ferroelectric  nature  of  the  materials,  they  also  exhibit  varying  degrees  of  hysteresis 
and  nonlinear  saturation  effects  at  moderate  to  high  drive  levels  as  illustrated  in  Figure  1.  As  detailed  in 
[19]  and  references  therein,  this  form  of  hysteresis  is  generally  attributed  to  the  impediment  of  domain  wall 
movement  by  material  inclusions  and  stress  nonhomogeneities  inherent  to  the  materials.  In  the  absence 
of  an  applied  field,  domain  walls  form  at  these  pinning  sites  to  minimize  the  associated  potential  energy. 
Various  experimental  investigations  have  illustrated  that  at  low  input  field  levels,  domain  wall  movement 

•Center  for  Research  in  Scientific  Computation,  North  Carolina  State  University,  Raleigh,  NC  27695  (rsmith@eos.ncsu.edu). 
This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  Number  NAS1- 
97046  while  this  author  was  a  consultant  at  the  Institute  for  Computer  Applications  in  Science  and  Engineering  (ICASE),  M/S 
132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 

t  ICASE,  M/S  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199  (z.ounaies@larc.nasa.gov).  This  research 
was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  Number  NAS1-97046  while  this 
author  was  in  residence  at  ICASE,  M/S  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199. 


1 


0.5 1 


_o.sU - . - i - 1 - « - ‘ - 

-6  -4  -2  0  2  4  6 

Electric  Field  (MV/m)  , 

Fig.  1.  Hysteresis  measured  in  a  PZT5A  wafer  in  response  to  a  1600  V  input. 

is  reversible  and  can,  at  least  conceptually,  be  attributed  to  the  bending  of  domain  walls  [3,  5,  14].  For 
higher  input  fields,  the  local  energy  barriers  associated  with  pinning  sites  are  overcome  and  domain  walls 
move  for  extended  distances  [3].  This  translation  of  domain  walls  across  pinning  sites  provides  an  irreversible 
mechanism  contributing  to  the  hysteresis  observed  in  ferroelectric  materials. 

The  hysteresis  inherent  to  ferroelectric  materials  can  be  accommodated  through  a  variety  of  techniques. 
The  simplest  means  for  minimizing  hysteresis  is  to  restrict  the  input  fields  or  stresses  to  sufficiently  low  levels 
to  maintain  quasilinear  behavior.  However,  this  severely  limits  the  capabilities  of  the  materials  and  is  not 
feasible  in  many  high  performance  applications  where  the  materials  are  advantageous.  For  certain  control  and 
damping  applications,  the  deleterious  effects  of  hysteresis  can  be  minimized  either  indirectly  through  feedback 
mechanisms  which  maintain  the  system  at  low  input  levels  or  directly  through  techniques  such  as  feedback 
linearization.  While  successful  in  certain  regimes,  both  approaches  are  often  significantly  handicapped  by 
the  phase  lags  associated  with  the  hysteresis  loops  [8, 17].  Furthermore,  open  loop  applications  which  require 
a  high  degree  of  accuracy  (e.g.,  micropositioning)  are  not  typically  amenable  to  these  feedback  techniques, 
and  hysteresis  must  be  accommodated  through  additional  design  criteria  or  models  which  can  be  used  to 
compensate  for  the  hysteresis  and  saturation  nonlinearities.  r 

In  this  paper,  we  consider  a  hysteresis  model  for  piezoelectric  materials  which  is  based  on  the  quantifi¬ 
cation  of  domain  and  domain  wall  mechanisms  inherent  to  the  materials.  The  model  is  based  on  the  theory 
developed  in  [18,  19]  for  the  general  characterization  of  hysteresis  in  ferroelectric  materials.  We  focus  here  on 
the  modeling  of  hysteresis  in  the  relation  between  the  applied  field  and  resulting  polarization;  the  resulting 
strains  can  then  be  specified  through  linear  constitutive  relations  in  the  manner  described  in  [18,  19].  In  the 
first  step  of  the  model  development,  the  hysteresis-free,  or  anhysteretic,  relation  between  an  applied  field 
and  the  resulting  polarization  is  quantified  through  three  techniques.  The  first  two  anhysteretic  models  are 
empirical  in  nature  and  are  used  to  provide  initial  estimates  for  parameters  in  the  final  hysteresis  model. 
The  third  is  based  on  the  classical  application  of  Boltzmann  statistics  and  provides  constitutive  relations 
which  are  applicable  at  low  drive  levels.  Hysteresis  is  then  modeled  through  the  characterization  of  the  re¬ 
versible  and  irreversible  motion  of  domain  walls  pinned  at  inclusions  in  the  material.  The  combination  of  the 
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components  provides  an  ODE  model  which  incorporates  the  nonlinear  constitutive  relations  and  hysteresis 
observed  in  various  PZT  compounds  at  high  drive  levels. 

While  the  model  is  based  on  the  theory  developed  in  [18, 19],  the  contributions  of  this  paper  are  threefold, 
(i)  The  first  is  the  extension  of  the  theory  to  include  nonsymmetric  hysteresis  loops.  The  resulting  model 
contains  the  symmetric  model  developed  in  [18,  19]  as  a  special  case  but  accommodates  the  nonsymmetric 
effects  due  to  poling  in  hard  PZT  materials,  (ii)  The  second  contribution  is  an  extensive  validation  of 
the  theory  for  commonly  employed  piezoelectric  materials  under  a  variety  of  drive  levels.  The  model  in 
[18, 19]  was  illustrated  in  the  context  of  the  relaxor  ferroelectric  PMN-PT-BT  employed  at  one  drive  level  at 
sufficiently  low  temperatures  for  it  to  be  ferroelectric.  In  this  paper,  the  predictive  capabilities  of  the  model 
are  illustrated  by  identifying  parameters  at  one  drive  level  and  then  using  the  resulting  model  to  predict 
the  PZT  material  behavior  at  other  input  levels.  The  model’s  capability  for  prediction  is  due  to  its  basis  in 
energy  principles  and  provides  it  with  an  important  advantage  in  broadband  applications.  Furthermore,  this 
capability  is  illustrated  for  the  commonly  employed  compounds  PZT5A,  PZT5H  and  PZT4.  (iii)  The  third 
contribution  of  the  paper  is  the  development  of  a  method  for  approximating  the  five  required  parameters 
in  the  model.  This  algorithm  is  analogous  to  that  developed  in  [13]  for  magnetic  materials  and  provides 
initial  estimates  for  the  parameters  through  a  comparison  of  the  model  with  physical  attributes  of  the  data 
including  the  coercive  field,  the  differential  susceptibility  at  various  points,  and  the  saturation  characteristics 
of  the  material.  These  estimates  can  be  employed  in  the  final  model,  if  sufficiently  accurate,  or  used  as  initial 
values  in  a  least  squares  fit  to  measured  data.  In  combination,  these  contributions  illustrate  the  applicability 
of  the  theory  for  piezoelectric  materials  and  extend  its  practical  feasibility  for  more  general  ferroelectric 
applications. 

We  note  that  the  current  model  is  quasistatic  and  isothermal  in  nature.  Moreover,  it  is  theoretically 
limited  to  materials  in  which  crystalline  anisotropies  are  not  significant.  While  initial  investigations  indicate 
that  its  applicability  extends  beyond  these  regimes,  such  applications  should  be  considered  with  caution 
until  the  underlying  physics  is  incorporated  in  the  models.  The  extensions  of  the  theory  to  accommodate 
frequency  and  thermal  effects  as  well  as  crystalline  anisotropies  are  under  investigation. 

A  brief  review  of  existing  models  for  hysteresis  in  piezoelectric  materials  is  summarized  in  the  remainder 
of  this  section  and  the  hysteresis  model  is  then  outlined  in  Section  2.  The  relationship  of  the  model  parameters 
to  physical  properties  exhibited  by  the  materials  is  detailed  in  Section  3  and  an  algorithm  for  estimating 
the  parameters  is  provided.  This  section  also  includes  a  discussion  of  least  squares  methods  which  can  be 
employed  for  final  parameter  determination.  In  Section  4,  the  model  is  fit  to  data  from  PZT5A,  PZT5H 
and  PZT4  wafers  under  a  variety  of  drive  conditions.  This  illustrates  both  the  accuracy  of  the  model  and 
its  capability  for  predicting  the  polarization  due  to  changing  drive  levels. 

1.1.  Existing  Models.  Hysteresis  models  for  piezoelectric  materials  can  be  roughly  categorized  as 
being  microscopic,  macroscopic  or  semi-macroscopic  in  nature.  Microscopic  theories  are  considered  as  those 
based  on  quantum  principles,  classical  elasticity  or  electromagnetic  relations,  or  thermodynamic  laws  applied 
at  the  lattice  or  grain  level.  While  such  theories  are  founded  on  the  underlying  physics,  they  often  require  a 
large  number  of  parameters  and  involve  states  which  are  difficult  or  impossible  to  measure  [15].  Moreover,  it 
is  difficult  to  incorporate  attributes  such  as  grain  boundaries  and  intergranular  stresses  in  the  models.  For 
this  reason,  microscopic  models  are  currently  limited  to  simple  stoichiometries  and  are  difficult  to  implement 
in  control  designs  due  to  the  large  number  of  required  parameters. 

Macroscopic  models  are  based  on  phenomenological  or  empirical  principles  and  are  considered  advanta¬ 
geous  when  the  underlying  physics  is  poorly  understood  or  difficult  to  characterize.  This  category  includes 
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models  based  on  shifted  anhysteretic  curves  which  provide  an  envelop  for  the  hysteresis  curves  (e.g.,  [23])  as 
well  as  various  Preisach  models  for  hysteresis  in  piezoelectric  materials  [6,  7,  8].  In  addition  to  the  generality 
provided  by  such  approaches,  the  resulting  models  can  for  certain  formulations  be  inverted  to  provide  com¬ 
pensators  for  linear  control  design  [6,  8,  21].  This  has  led  to  the  fairly  widespread  use  of  Preisach  techniques 
for  modeling  hysteresis  in  a  variety  of  smart  materials  including  piezo  ceramics.  The  disadvantage  with  this 
approach  lies  in  the  fact  that  it  is  difficult  to  employ  known  physics  or  physical  measurements  to  directly 
construct  the  model  or  update  the  model  parameters  to  accommodate  changing  operating  conditions.  Fur¬ 
thermore,  to  provide  flexibility  for  a  variety  of  drive  conditions  (e.g.,  broadband  or  transient  conditions), 
it  is  necessary  to  employ  a  large  number  of  nonphysical  parameters  which  makes  real-time  implementation 
difficult. 

Semi-macroscopic  theories  are  derived  using  a  combination  of  these  approaches.  They  typically  employ 
energy  relations  to  characterize  attributes  of  the  polarization  switching  mechanism  and  then  use  macroscopic 
averages  to  obtain  parameters  and  subsequent  models  for  the  bulk  behavior  of  the  polycrystalline  material. 
For  example,  the  theory  of  Chen  and  Lynch  [2]  employs  energy  relations  to  quantify  the  polarization  and 
strain  at  the  grain  level.  Macroscopic  averaging  over  the  grains  is  then  used  to  characterize  the  aggregate 
behavior  of  the  material.  A  model  for  quasistatic  hysteresis  was  obtained  by  Huang  and  Tiersten  [9]  through 
the  incorporation  of  internal  variables  in  the  thermodynamic  state  to  incorporate  the  irreversible  effects 
inherent  to  domain  switching.  The  resulting  model  employs  seven  nonlinear  material  coefficients  and  three 
saturation  polarization  coefficients  to  provide  a  characterization  of  the  hysteresis  in  the  relations  between 
the  field  and  polarization  and  field  and  strain. 

The  model  employed  in  this  paper  fits  in  the  latter  category.  Electrostatic  energy  relations  are  employed 
to  quantify  the  reversible  and  irreversible  effects  of  domain  wall  bending  and  translation.  Macroscopic  aver¬ 
ages  then  provide  a  model  whose  five  parameters  quantify  the  bulk  attributes  of  the  material.  Construction 
in  this  manner  provides  the  model  with  its  predictive  capabilities  and  facilitates  the  determination  and 
updating  of  parameters  from  measured  data.  Both  aspects  prove  advantageous  in  control  applications. 

This  provides  an  overview  of  certain  models  which  specifically  address  the  modeling  of  hysteresis  in 
piezoelectric  materials.  The  reader  is  referred  to  [18,  19]  for  additional  discussion  of  previous  hysteresis 
models  for  general  ferroelectric  materials. 

2.  Hysteresis  Model.  We  summarize  in  this  section  those  aspects  of  the  hysteresis  theory  presented 
in  [18,  19]  which  are  relevant  to  piezoelectric  materials.  In  the  first  step  of  this  development,  various  models 
for  the  anhysteretic  polarization  are  outlined.  Hysteresis  is  then  incorporated  through  the  quantification  of 
the  energy  required  to  bend  and  translate  domain  walls. 

2.1.  Anhysteretic  Polarization.  The  anhysteretic  polarization  curve  represents  the  minimal  energy 
states  that  would  be  attained  by  the  material  in  the  absence  of  inclusions  or  imperfections.  As  illustrated 
by  the  modeled  curves  in  Figure  2,  the  anhysteretic  has  a  burst  region  near  the  origin,  where  dipole  rotation 
toward  preferred  ionic  configurations  produces  large  changes  in  polarization,  and  exhibits  saturation  at  high 
field  levels  where  charge  distributions  prohibit  further  changes.  In  contrast  to  the  measured  polarization 
curves  at  high  drive  levels,  the  anhysteretic  polarization  is  single- valued  and  reversible. 

Three  techniques  are  outlined  for  modeling  the  anhysteretic  curve.  The  first  two  are  empirical  and 
provide  expressions  which  are  used  to  specify  parameters  in  the  algorithms  developed  in  Section  3.  The  last 
technique  specifies  the  probability  of  dipoles  occupying  certain  energy  states  through  classical  Boltzmann 
statistics.  The  resulting  representations  are  used  as  kernels  in  the  hysteresis  model. 


4 


-0.2  0  0.2 
Electric  Field  (MV/m) 


Fig.  2.  Ising  spin  and  Langevin  models  for  the  anhysteretic  polarization. 


The  first  model  is  based  on  the  extension  of  the  Frolich-Kennelly  anhysteretic  model  for  magnetic 
materials  [11,  page  94]  to  ferroelectric  materials.  This  yields  the  expression 

.  7  E 

an  l  +  7£ 

where  E  and  Pan  respectively  denote  the  applied  electric  field  and  anhysteretic  polarization.  One  parameter 
can  be  eliminated  by  enforcing  the  saturation  behavior  Pan  -»•  Ps  as  E  ->  oo,  where  P8  is  the  saturation 
polarization.  This  then  yields  the  relation 

/o  D  p  — 

C2*1)  i“n  1  +  7E  * 

While  not  employed  in  the  final  hysteresis  model,  the  expression  (2.1)  provides  a  means  of  estimating  certain 
parameters  in  the  model. 

A  second  model  for  the  anhysteretic  polarization  is  the  empirical  expression 
(221  P  - 

<M)  ’  ./T+W‘ 

developed  by  Piquette  and  Forsythe  [16].  The  initial  behavior  and  slope  of  this  curve  at  high  fields  is  similar 
to  that  of  (2.1)  while  the  mid-range  behaviors  and  saturation  values  differ.  The  expression  (2.2)  will  also  be 
used  to  specify  initial  parameter  values  in  the  algorithms. 

The  third  model  employs  Boltzmann  statistics  to  specify  the  probability  of  dipoles  occupying  certain 
energy  states.  As  detailed  in  [18,  19],  the  balance  of  thermal  and  electrostatic  energies  while  employing 
the  assumption  that  the  material  is  isotropic  and  the  orientation  of  cells  can  be  in  any  direction  yields  the 
Langevin  equation 


=  P8  coth 


(Ee\  _  a_' 

\  a  )  Ee ' 


Ee  =  E  +  aPa , 


for  the  anhysteretic  polarization.  Here 
(2.4) 
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denotes  the  effective  field  acting  on  dipole  moments  in  the  material  while  a  =  E/Pa  where  E  denotes  a 
scaling  electric  field.  The  constant  a  is  specified  by  a  =  where  T  and  Tc  are  the  operating  and  Curie 
temperatures  for  the  material.  Because  the  scaling  field  E  is  unknown,  the  constants  a  and  a  are  estimated 
using  either  the  algorithm  derived  in  Section  3  or  a  least  squares  fit  to  data. 

A  second  expression  derived  using  this  Boltzmann  approach  is  the  Ising  spin  relation 


(2.5) 


=  P’tanh(f) 


which  results  from  the  assumption  that  the  dipole  moments  can  occupy  only  two  discrete  orientations:  in  the 
direction  of  the  applied  field  or  opposite  to  it.  As  noted  in  Figure  2  where  the  relations  (2.3)  and  (2.5)  are 
compared,  the  restrictions  on  possible  moment  orientations  causes  the  Ising  spin  relation  to  reach  a  higher 
saturation  state  than  observed  for  the  Langevin  model  with  equivalent  inputs.  Further  discussion  concerning 
the  merits  of  the  two  relations  as  well  as  details  regarding  their  derivation  are  provided  in  [18,  19]. 

The  relations  (2.3)  and  (2.5),  with  the  effective  field  specified  by  (2.4),  yield  anhysteretic  curves  which 
are  rotationally  symmetric  about  the  origin  as  illustrated  in  Figure  2.  For  PMN-PT-BT  or  soft  piezoceramic 
materials,  these  expressions  are  sufficient  since  the  relation  between  the  applied  field  and  resulting  polar¬ 
ization  is  also  rotationally  symmetric.  This  is  not  the  case  for  hard  PZT  materials  which  have  been  poled, 
however,  and  the  expressions  must  be  modified  to  incorporate  the  observed  asymmetries. 

To  motivate  modifications  which  incorporate  the  observed  asymmetries,  consider  the  data  from  a  poled 
PZT4  wafer  which  is  plotted  in  Figure  3.  It  is  observed  that  the  poling  produces  a  bias  polarization  P0  and 
field  E0  which  causes  the  positive  and  negative  remanence  points  and  coercive  fields  to  differ.  These  values 
bias  the  effective  field  so  we  consider  the  modified  expression 


(2.6) 


Ee  =  (E-Eo)  +  a(P-P0). 


In  this  case,  the  effective  field  reflects  the  bias  field  in  the  absence  of  an  applied  field  E. 

To  accommodate  both  asymmetric  major  loops  and  biased  minor  loops,  we  also  consider  the  modified 
anhysteretic  expressions 

(2.7)  Pan  =  Pl  +  Ps  [coth 
and 

(2.8)  Pan  =  Pi  +  P8tanh 

where  Pi  is  a  scaling  polarization.  The  curves  produced  by  the  Langevin  model  (2.7)  with  Pi  =  0  and 
Pl  =  pQ  ^  0  are  illustrated  in  Figure  4.  In  both  cases,  the  curves  are  asymmetric  with  respect  to  the  origin. 
For  Pi  =  0,  the  curve  is  rotationally  symmetric  about  a  point  ( E ,  P)  determined  through  the  numerical 
solution  of  (2.7)  or  (2.8)  whereas  for  Pi  -  Po,  the  change  of  variables  P  =  P-Po,£  =  P-Po  reveals  that 
the  curve  is  rotationally  symmetric  about  the  bias  point  (Eq^Pq). 
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Fig.  3.  Data  from  a  poled  PZT4  wafer  in  response  to  a  2000  V  input 


Electric  Field  (MV/m)  Electric  Field  (MV/m) 

(a)  (b) 

Fig.  4.  Anhysteretic  curves  produced  by  the  Langevin  model  (2.7)  with  a  =  3.4  x  106  Vm/C  and 
a  =  8.0  x  105  C/m2;  (a)  E0  =  -0.5  x  106  V/m,  P0  =  .05  C/m 2  and  Pi  =  0,  (b)  E0  =  -0.5  x  106  V/m, 
P0  =  .05  C/m 2  and  Px  =  P0. 

2.2.  Domain  Wall  Model.  The  relations  (2.7)  or  (2.8)  can  be  used  to  model  the  observed  polarization 
at  low  drive  levels  but  are  inappropriate  at  moderate  to  high  drive  levels  since  they  do  not  incorporate  the 
hysteresis  inherent  to  the  materials.  As  detailed  in  [19]  and  included  references,  sigmoidal  hysteresis  of  the 
type  depicted  in  Figure  1  is  typically  attributed  to  the  energy  required  to  translate  domain  walls  across 
pinning  sites  in  the  material.  At  low  field  levels,  the  walls  remain  close  to  the  equilibrium  position  and 
the  motion  is  reversible.  From  an  energy  perspective,  the  variations  are  not  sufficient  to  cross  a  barrier  in 
the  potential  well.  The  motion  becomes  irreversible  when  sufficient  energy  is  provided  to  cross  the  potential 
barrier.  Physically,  this  can  occur  when  the  domain  wall  intersects  a  remote  pinning  site  and  is  the  mechanism 
underlying  domain  wall  translations.  The  resulting  irreversible  polarization  PiTT  and  reversible  polarization 


7 


Prev  are  then  summed  to  obtain  the  total  polarization.  This  approach  follows  that  employed  by  Jiles  and 
Atherton  in  their  corresponding  hysteresis  model  for  ferromagnetic  materials  [12]. 

To  quantify  the  irreversible  polarization,  it  is  noted  in  [18,  19]  that  the  polarization  level  for  a  given 
effective  field  can  be  expressed  as  that  for  the  ideal  case  minus  losses  required  to  break  pinning  sites.  This 
yields  the  relation 

_  ,  dPirr 

(2.9)  Pirr  —  Pan  ~~  &  • 

The  parameter  k  is  defined  by  k  =  where  n  denotes  the  average  density  of  pinning  sites,  {£„)  is 

the  average  energy  for  180°  walls  and  p  is  an  average  dipole  moment.  Because  the  density  and  energy  of 
individual  pinning  sites  are  unknown,  the  parameter  k  must  be  estimated  for  a  given  material. 

The  formulation  of  (2.9)  in  terms  of  the  applied  field  E  yields  the  differential  equation 


dPirr  _  Pan  Pjrr _ 

dE  Sk  Oi{Pa  ji  Pirr ) 


specifying  the  irreversible  polarization.  The  parameter  S  =  sign(d-E)  ensures  that  the  energy  required  to 
break  pinning  sites  always  opposes  changes  in  polarization.  As  discussed  in  [18,  19],  while  this  expression 
is  adequate  in  most  regimes,  it  can  yield  nonphysical  solutions  when  the  field  is  reversed  from  saturation 
for  materials  which  exhibit  significant  hysteresis  and  are  driven  at  high  levels.  The  enforcement  of  solely 
reversible  polarization  changes  in  this  regime  eliminates  this  discrepancy  and  yields  the  relation 


(2.10) 


dPirr  ^ _ P an  Pjrr _ 

dE  kS  ~  Ct  {Pan  “T  Pirr ) 


where 

~  (  1  ,  {dE  >  0  and  P  <  Pan }  or  {dE  <  0  and  P  >  Pan } 

[  0  ,  otherwise . 

The  second  component  of  the  polarization  is  the  reversible  polarization  which  models  the  effects  of 
domain  wall  bending.  To  first  approximation,  this  is  modeled  by  the  relation 


(2.11)  Prev  —  c(-P 2  n  Pirr ) 

where  c  is  a  parameter  which  must  be  estimated  for  the  specific  application  (see  [18,  19]). 
The  total  polarization  is  then  given  by 


P  —  Prev  +  Pirr 

or  equivalently, 

(2.12)  P  =  cPan  +  (1  -  c)Pirr  . 

To  implement  the  model,  the  effective  field  for  a  given  field  and  irreversible  polarization  level  is  computed 
using  (2.6).  This  effective  field  value  is  then  employed  in  either  (2.7)  or  (2.8)  to  compute  the  corresponding 
anhysteretic  polarization.  The  subsequent  irreversible  polarization  is  determined  by  numerically  integrating 
(2.10).  The  total  polarization  is  then  specified  by  (2.12).  The  determination  of  parameters  is  addressed  in 
the  next  section  while  the  prediction  capabilities  of  the  model  are  illustrated  in  Section  4. 
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3.  Determination  of  Parameters.  The  implementation  of  the  model  requires  the  determination  of 
the  parameters  a,  a,  k ,  c  and  Ps  for  a  given  material.  Recall  that  a  quantifies  the  amount  of  dipole  coupling  in 
the  effective  field  while  a  incorporates  the  relative  thermal  effects  which  are  balanced  with  the  electrostatic 
energy  to  model  the  anhysteretic  polarization.  As  detailed  in  [19],  increasing  a  or  decreasing  a  leads  to 
increased  slopes  in  the  anhysteretic  curve  and  corresponding  polarization.  The  parameter  k  is  a  macroscopic 
average  of  the  energy  required  to  break  pinning  sites.  Hence  large  k  values  are  associated  with  wider 
hysteresis  loops  (e.g.,  hard  PZT  materials).  The  parameter  c  quantifies  the  average  degree  to  which  domain 
walls  bend  before  translating  across  pinning  sites;  hence  it  too  will  be  larger  for  hard  materials  than  for 
soft.  Finally,  Ps  denotes  the  theoretical  saturation  value  beyond  which,  polar  interactions  prevent  further 
increases  in  polarization.  This  definition  is  in  accordance  with  the  corresponding  definition  for  the  saturation 
magnetization  (e.g.,  see  [11])  and  should  not  be  confused  with  the  definition  employed  in  several  texts  (e.g., 
[10,  page  36])  for  the  saturation  value  obtained  by  extending  the  slope  at  tip  reversal  back  through  the 
vertical  axis.  Further  details  regarding  the  derivation  and  physical  interpretation  of  these  parameters  are 
provide  in  [18]. 

In  this  section,  we  present  two  methods  for  determining  the  parameters  based  on  data  measurements 
from  a  given  material.  The  first  uses  the  measured  values  of  the  anhysteretic  and  initial  susceptibilities 
(if  available),  the  measured  remanence  polarization,  the  coercive  field  and  peak  tip  polarization,  as  well 
as  the  differential  susceptibilities  at  these  points,  to  provide  constraints  which  permit  the  determination  of 
the  parameters.  This  approach  is  analogous  to  that  employed  in  [13]  for  magnetic  materials  but  leads  to 
a  different  algorithm  for  determining  the  parameters.  This  technique  highlights  the  physical  nature  of  the 
parameters  and  is  direct  to  implement,  but  can  lead  to  model  fits  with  limited  accuracy  since  it  employs 
minimal  information  concerning  the  hysteresis  curve.  The  second  approach  determines  the  parameters 
through  a  least  squares  fit  to  the  data.  This  provides  highly  accurate  model  fits,  but  requires  fairly  good 
initial  guesses  for  the  parameters  to  reach  optimal  values.  In  practice,  we  employ  the  first  technique  to 
obtain  initial  parameter  values.  In  many  cases,  the  resulting  model  fit  is  satisfactory.  If  refinement  is 
necessary,  however,  these  parameter  values  can  be  employed  as  initial  estimates  in  the  least  squares  routine. 
In  combination,  the  two  approaches  provide  a  systematic  and  robust  means  of  determining  the  parameters 
for  a  given  material. 

3.1.  Direct  Parameter  Determination  from  Experimental  Data.  The  saturation  polarization 
Ps  has  the  most  direct  physical  connotation  and  can  be  estimated  from  the  data  at  high  (near  saturation) 
drive  levels.  These  estimates  can  then  be  refined  using  either  of  the  following  methods. 


3.1.1.  Differential  Susceptibility  Relations.  To  determine  the  constraints  used  to  specify  the  re¬ 
maining  parameters,  we  consider  the  differential  susceptibility  at  various  points  in  the  hysteresis  cycle.  From 
the  definition  (2.10)  for  the  differential  susceptibility  of  the  irreversible  polarization  and  (2.12)  for  the  total 
polarization,  we  see  that  two  cases  need  to  be  considered  when  computing  When  the  field  is  first  reversed 
from  saturation,  the  only  changes  in  polarization  are  due  to  the  reversible  effects  of  domain  wall  relaxation. 
This  motivates  the  inclusion  of  the  switch  S  and  yields  the  expression 


(3-1) 


dP  _  dPan 
dE~C  dE 


for  the  differential  susceptibility  in  that  region.  For  the  remainder  of  the  hysteresis  cycle,  the  combination 
of  (2.10)  and  (2.12)  yields  the  expression 


dP  Pan  Pirr  ,  dP 2n 

1e~K  ~C)k5-  a(Pan  -  Pirr)  C~dE~ 
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A  second  application  of  the  polarization  relation  (2.12)  can  be  used  to  eliminate  the  irreversible  polarization 
which  yields 

(o  o')  — -(1-c) _ Pan  -  P _ +  cdPan 

(3-2)  dE  ~  {  ’  kS(  1  -  c)  -  a(Pan  -P)  dE  • 

The  evaluation  of  (3.2)  requires  the  determination  of  which  depends  upon  the  model  being  em¬ 
ployed.  For  the  examples  in  Section  4,  we  employ  the  Langevin  expressions  (2.3)  or  (2.7)  and  we  illustrate 
with  (2.3)  here.  For  general  polarization  values,  the  effective  field  is  Ee  =  E  +  aP  which  yields 


^  =  +  -csch2 

dE  a  V  dE  I 


E  +  aP 


For  the  single- valued  global  anhysteretic  curve  depicted  in  Figure  2,  the  effective  field  is  specified  by  (2.4) 
which  yields  the  implicit  relation 

^  fl  +  a^-)  [-csch2  (- -  +  ( — — —5—)  • 

dE  a  \  dE  J  [  \  a  )  \E  +  aPanJ. 

Initial  Susceptibilities 

Taking  the  limits  as  E  ->  0,  Pan  ->  0,  and  letting  Xan  denote  the  differential  susceptibility  at  the  origin, 
as  depicted  in  Figure  5,  yields 

Ps 

Xan  =  (1  +  aXan) 


For  given  Ps,  this  yields  the  expression 


3a  -  aPs 


3  aXon  P i s 


relating  a  to  a. 

The  second  characteristic  which  can  be  employed  at  the  origin  is  the  initial  differential  susceptibility 
Xm-  Taking  the  limits  E  0,  P  0  in  (3.2)  yields  the  relation 

(3.6)  Xin  ~  cXon 

which  can  be  used  to  specify  c  if  x%n  and  Xan  can  be  measured  or  approximated.  To  evaluate  the  slope  of 
the  initial  polarization  curve,  we  take  the  limit  E  -»  0,P  -*  0  in  (3.2)  with  the  anhysteretic  specified  by 
(3.3).  This  yields 

P  c 

X in  “  H"  axin) 


from  which  it  follows  that 


3aXin  cP 8 

cEsXin 


Equating  (3.5)  and  (3.7)  and  employing  (3.6)  to  eliminate  Xm  then  yields  the  expression 
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Fig.  5.  Hysteresis  curve  with  differential  susceptibilities  employed  for  parameter  determination . 

Once  Xan  has  been  approximated  or  determined  from  experimental  data,  (3.8)  can  be  used  to  approximate 
the  parameter  a  and  a  can  be  determined  from  (3.5). 


Susceptibility  at  Field  Reversal 


We  consider  next  the  behavior  of  the  hysteresis  loop  near  the  tip  value  (Emi  Pm).  It  was  previously  noted 
that  directly  after  field  reversal,  the  differential  susceptibility  is  dependent  only  on  domain  wall  relaxation 
which  yields  the  expression  (3.1).  Furthermore,  if  the  polarization  is  sufficiently  close  to  saturation,  the 
differential  susceptibility  can  be  approximated  by 


(3.9) 


dP  _  dPan 
dE~~dE~ 


in  the  region  before  field  reversal  (see  also  [13]  for  the  magnetic  case).  Letting  xti  and  X™  respectively  denote 
the  differential  susceptibilities  before  and  after  field  reversal  at  the  tip  loop  (see  Figure  5),  the  combination 
of  (3.1)  and  (3.9)  yields  the  expression 


(3.10) 


for  the  reversibility  coefficient.  We  note  that  the  accuracy  of  this  expression  improves  for  measurements  near 
saturation. 

Furthermore,  since  S  =  1  before  field  reversal,  the  consideration  of  the  approximation  (3.9)  in  (3.2)  yields 


(3.11) 


X' 


+  _ 


P an  {.Era)  Pm 


&5(1  —  c)  —  Oi(Pan(Em)  P m) 


where,  for  the  symmetric  case,  Pan(Em)  is  specified  by  (2.3)  or  (2.5)  with  Ee  =  Em  +  aPm.  The  expression 
(3.11)  can  be  solved  explicitly  for  k  or  implicitly  for  a  or  a. 


ll 


Remanence  Point 


Similar  constraints  are  provided  by  matching  the  differential  susceptibilities  at  the  coercive  and  rema¬ 
nence  points.  At  remanence,  <5  =  —  1,  E  =  0  and  P  =  Pr  so  that  the  relation  (3.2)  has  the  form 

(312)  -V  =  (1  —  c) _ -Pqn(fi)  -  Py -  dPan{Pr) 

(3‘12)  Xr  {  ’  k5(  1  -  c)  -  a(Pon(Pr)  -  Pr)  dE 

where,  for  the  symmetric  Langevin  expression  (2.3), 


ftn(ft)  =  ft  COth 


(  «ft\ _ 

\  a  )  aPr 


dP an 
dE 


(ft)  = 


—  (1  +  ax  r) 
a 


Because  the  relation  (3.12)  is  implicit  in  the  variables  a  and  a,  it  is  solved  through  either  root-finding  or 
minimization  techniques. 


Coercive  Field 

At  the  coercive  field,  S  =  1,  E  =  Ec  and  P  =  0,  so  the  differential  susceptibility  is 

,,  '  Panm  .  „dPan(Ec) 

Xc  -  U  C) kS(  1  -  c)  -  aPan(Ec)  +  dE  • 

This  can  be  solved  for  k  to  yield 
(3.13)  k  =  Pan{Ec)  YZTC  + 

For  specified  values  of  a,a,c  and  P6 ,  (3.13)  provides  an  estimate  for  the  average  energy  required  to  break 
pinning  sites  in  both  hard  and  soft  materials. 

As  detailed  in  [11,  pages  170-171]  for  magnetic  materials,  this  expression  can  be  simplied  significantly 
for  soft  materials  when  the  reversible  component  is  negligible  and  hence  the  approximation  c  =  0  is  valid. 
For  materials  with  low  coercivity,  it  has  also  been  observed  that  the  differential  susceptibility  at  the  coercive 
point  is  approximately  equal  to  the  slope  of  the  anhysteretic  curve  at  the  origin  so  that 


As  will  be  illustrated  in  the  examples,  the  values  predicted  by  (3.15)  can  be  used  as  initial  values  for  the 
parameter  estimation  routines  and  in  many  cases  are  quite  accurate. 
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3.1.2.  Approximation  of  Xan  and  Xin ■  The  relations  (3.4)  -  (3.13)  provide  constraints  which  will 
be  used  for  constructing  an  algorithm  for  determining  the  parameters  a,  a,  c and  k.  For  cases  in  which  the 
measured  values  of  Xan  and  \in  are  available,  the  constraints  (3.4)  and  (3.6)  can  be  applied  directly.  In 
applications  in  which  the  material  starts  from  a  poled  state,  however,  these  values  may  be  unavailable  which 
necessitates  their  approximation. 

For  ferromagnetic  materials  whose  hysteresis  loops  are  rotationally  symmetric  about  the  origin,  it  has 
been  observed  that  for  a  given  field  value,  points  Min  on  the  initial  magnetization  curve  can  be  approximated 
by  the  average 

Min  =  (M1  +  M2)/2 


where  Mi  and  M2  are  the  corresponding  magnetizations  on  the  upper  and  lower  hysteresis  curves  [1,  page 
511],  [4].  A  corresponding  assumption  for  ferroelectric  materials  yields  the  approximation 

(3.16)  Zm  =  (Xr  +  Xr)/2 


where  Xr  and  x7  denote  the  differential  susceptibilities  at  the  positive  and  negative  remanence  points  (see 
Figure  5). 

The  initial  slope  of  the  anhysteretic  curve  can  be  approximated  using  the  alternative  model  expressions 
(2.1)  and  (2.2).  To  eliminate  the  single  parameter  appearing  in  each  expression,  we  enforce  the  condition 


dPgn 

~dE 


{Pm) —  Xm 


where  =  (x+  +  Xm)/2*  In  expression  (2.1),  this  yields 

—  ~(2-^mXm  —  Ps)  y/Ps  ~~  ^ sPmXm 
^  ^XmPm 


and  the  corresponding  initial  slope 


Xani  =  7 Ps  • 


The  corresponding  initial  slope  predicted  by  (2.2)  is 

Xari2  =  y/ftPs 


where 

_  1  61/3  2  P2S  1 

P~  GElxm  +  E*mXmb'l*  El 

with 

b  =  -12 P2S  [9 EmXm  -  vV-4 p!  +  27Elxl]  . 

The  initial  slope  of  the  anhysteretic  is  then  approximated  by 
(3.17)  Xan  =  (Xani  XanT)!^' 

We  note  that  the  accuracy  of  the  approximations  (3.16)  and  (3.17)  is  dependent  upon  the  specific  material 
and  operating  conditions  under  consideration.  To  accommodate  inaccuracies  due  to  either  these  approxima¬ 
tions  or  measurement  errors,  the  following  algorithms  employ  the  remaining  constraints  to  refine  estimated 
parameter  values. 
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3.1.3.  Determination  of  Parameters.  The  susceptibility  and  saturation  criteria  can  be  combined 
to  provide  algorithms  for  specifying  the  parameters  a,  a,k  and  c  and  updating  P3.  In  certain  cases,  the 
relations  provide  implicit  constraints  on  the  variables  which  necessitates  the  use  of  root  finding  techniques 
(e.g.,  the  Matlab  routine  fzero.m)  for  solution.  Algorithms  1  and  2  differ  in  the  manner  through  which 
parameter  estimates  are  refined.  In  both  cases,  initial  estimates  are  obtained  directly  from  the  constraints. 
In  Algorithm  1,  these  estimates  are  refined  by  iterating  through  the  constraints.  While  direct  to  implement, 
this  technique  does  not  enforce  criteria  which  guarantee  convergence.  The  second  algorithm  is  more  robust 
since  refinement  is  accomplished  through  the  minimization  of  a  functional  which  simultaneously  incorporates 
all  constraints. 

The  accuracy  of  the  parameters  obtained  using  either  method  is  dependent  upon  the  degree  to  which 
the  slope  information  at  the  initial,  remanence,  coercive  and  extreme  points  quantifies  the  overall  behavior 
of  the  hysteresis  curve.  For  cases  in  which  this  information  is  not  sufficient,  the  parameter  estimates  can  be 
employed  as  initial  values  in  the  least  squares  routine  described  in  Section  3.2.  This  incorporates  the  full 
behavior  of  the  hysteresis  curve  and  produces  model  fits  which  are  optimal  in  a  least  squares  sense. 

Algorithm  1  (Iterative  Refinement): 

(A)  Determine  Initial  Parameter  Values: 

(1)  Specify  c :  From  (3.10), 

c=*f. 

Xm 

This  can  be  updated  using  c  -  Xin/Xan  from  (3.6)  using  either  measured  values  of  Xin  and  Xan 

or  approximates  given  by  (3.16)  and  (3.17),  respectively. 

(2)  Specify  a:  From  (3.8), 

Ps 

a  =  - - 

3 Xan 

where  Xan  is  either  measured  or  approximated  using  (3.17). 

(3)  Specify  a:  Solve  (3.12) 

\ _ PgnjPr)  ~~  Py _ .  dPan(Pr) 

Xr~(  C)kS(l-c)-a(Pan(PT)-Pr)+  dE 

using  k  —  Ec  from  (3.15). 

(4)  Specify  k :  From  (3.13) 

*  =  P“"(E‘>  +  ' 

(B)  Iterative  Refinement:  Iterate  until  convergence  is  achieved. 

(i)  Solve  (3.12)  for  a 

(ii)  Solve  (3.11)  for  a 

(iii)  Solve  (3.13)  for  k 
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Algorithm  2  (Simultaneous  Refinement): 

(A)  Determine  Initial  Parameter  Values:  Same  as  Algorithm  1. 

(B)  Simultaneous  Refinement:  Solve 


min||^)|| 


where  q  —  [a,a,c,k\  and  P(q)  =  [Fi(q),F2(g),Fs(q)]T  with 


Fl{q)  =  Xr-(l-c) 
E2  (q)  =  k  —  Pan(Fc) 

F3  (q)  =  x+ - 


Pan(Pr)  ~  Pr 


k8(l  -  C)  -  a(Pan(Pr)  ~  Pr) 
1 


—  C 


dPqnjPr) 
dE 


1  —  c 

Pan{Em)  ~~  Pi 


dPgn(Bc) 

Xc  c  dE 


k5(  1  “  C)  -  a(Pan{Em)  ~  Pm) 

Note:  One  can  also  consider  q  =  [a,  a,  c,  k,  P8]  to  additionally  update  Ps . 


(From  (3.12)) 
(From  (3.13)) 
(From  (3.11)) 


3.2.  Least  Squares  Determination  of  Parameters.  The  algorithms  developed  in  Section  3.1  high¬ 
light  the  physical  nature  of  the  parameters  and  are  straightforward  to  implement.  Because  they  incorporate 
a  limited  amount  of  information  through  the  constraints,  however,  they  may  not  provide  sufficient  accuracy 
in  certain  applications.  A  least  squares  method  of  the  type  described  here  incorporates  the  polarization 
values  measured  throughout  the  hysteresis  cycle  and  can  be  used  to  obtain  parameters  that  yield  models 
which  optimally  fit  the  data  in  a  least  squares  sense. 

To  formulate  the  least  squares  parameter  estimation  problem,  let  ( E{ ,  Pi),  i  —  1,  •  •  • ,  A/*,  denote  the  field 
and  corresponding  polarization  values  measured  throughout  the  hysteresis  cycle.  Furthermore,  let  P{Ei\q) 
denote  the  parameter-dependent  model  solutions  specified  by  (2.12).  For  admissible  parameters  q  €  Q,  we 
then  solve  the  optimization  problem 

A  ^,2 

(3.18)  min 2^  \P{Ei;q)  -  Pi\  . 

q  i=l 

An  initial  value  qQ  can  be  specified  either  through  a  priori  information  or  the  parameter  estimates  obtained 
using  the  algorithm  developed  in  the  previous  section.  Model  fits  obtained  using  this  procedure  are  provided 
in  the  next  section. 

4.  Model  Validation.  To  illustrate  the  performance  and  prediction  capabilities  of  the  model  and  pa¬ 
rameter  estimation  methods,  we  consider  the  characterization  of  hysteresis  in  a  variety  of  PZT  compounds. 
Specifically,  we  consider  its  performance  for  PZT5A,  PZT5H,  and  PZT4  wafers  as  well  as  its  flexibility  for 
characterizing  the  hysteresis  in  patches  having  different  geometries.  The  materials  and  geometric  config¬ 
urations  which  we  consider  are  summarized  in  Table  1.  All  data  was  collected  at  200  mHz  to  maintain 
quasistatic  operating  conditions  with  the  exception  of  a  1  Hz  set  for  PZT5A  which  is  included  to  illustrate 
that  even  at  1  Hz,  frequency-dependent  effects  are  observed  in  the  material. 

The  parameters  predicted  by  Algorithm  2  in  Section  3.1  and  the  least  squares  method  from  Section  3.2  are 
compiled  in  Table  2.  The  measured  field,  polarization  and  slope  characteristics  employed  in  the  algorithms 
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are  summarized  in  Table  3.  We  note  that  measured  values  of  the  initial  anhysteretic  and  normal  polarization 
curves  were  not  available  so  Xin  and  Xan  were  approximated  using  (3.16)  and  (3.17).  A  comparison  of  the 
parameters  predicted  by  Algorithm  2  and  the  least  squares  method  illustrates  that  while  they  are  close,  the 
least  squares  fit  refines  the  values  to  provide  better  fits  throughout  the  hysteresis  cycle. 

Example  1  (PZT5A): 

We  first  consider  the  characterization  of  hysteresis  exhibited  by  PZT5A  actuators  with  three  sets  of  data 
being  considered.  The  first  two  illustrate  that  different  actuators  can  exhibit  slightly  different  hysteresis 
characteristics  which  necessitates  the  updating  of  parameters.  The  third  illustrates  that  frequency-dependent 
effects  can  be  present  even  at  1  Hz  which  motivates  the  use  of  200  mHz  data  for  quasistatic  material 
characterization. 


Table  1.  Materials ,  frequencies  and  geometrical  configurations  considered  in  the  examples. 


Material 

Frequency 

Geometry 

Dimensions 

PZT5A 

200  mHz 

Circular 

2.54  cm  Diameter,  0.0254  cm  Thick 

Example  1 

PZT5A 

PZT5A 

200  mHz 

1  Hz 

Rectangular 

Rectangular 

1.7  cm  x  0.635  cm  x  0.0381  cm 

1.7  cm  x  0.635  cm  x  0.0381  cm 

Example  2 

PZT5H 

200  mHz 

Rectangular 

3.81  cm  x  0.635  cm  x  0.0381  cm 

Example  3 

PZT4 

200  mHz 

Rectangular 

3.81  cm  x  0.635  cm  x  0.0381  cm 

Table  2.  Parameters  determined  using  Algorithm  2  and  the  least  squares  method  of  Section  3.2  from 
data  collected  at  200  mHz.  PZT5A*  is  circular  and  PZT5A t  is  rectangular. 


PZT5A* 

Alg.  2 

(1600  V) 
Least  Sq. 

PZT5At  (1600  V) 
Alg.  2  Least  Sq. 

PZT5H  (2200  V) 
Alg.  2  Least  Sq. 

PZT4  (1800  V) 

Alg.  2  Least  Sq. 

a 

3.6  x  106 

3.6  x  106 

3.1  x  106 

3.7  x  106 

4.0  x  106 

4.2  x  10® 

6.5  x  10® 

6.4  x  10® 

a 

4.4  x  105 

4.2  x  10® 

4.2  x  10® 

4.1  x  105 

5.8  x  105 

6.4  x  10s 

8.3  x  10® 

8.0  x  106 

k 

1.9  x  10® 

1.8  x  106 

1.8  x  10® 

1.5  x  106 

1.1  x  106 

1.0  x  10® 

2.5  x  10® 

1.5  x  10® 

c 

0.18 

0.30 

0.22 

0.15 

0.14 

0.20 

0.37 

0.40 

Ps 

0.49 

0.49 

0.49 

0.49 

0.425 

0.425 

0.44 

0.44 

Actuator  1  (Circular  Patch): 

As  summarized  in  Table  1,  the  circular  wafer  had  a  diameter  of  2.54  cm  and  was  0.0254  cm  (10  mils) 
thick.  The  wafer  was  depoled  before  use  so  the  initial  cycle  contained  transient  behavior  as  the  material  was 
polarized.  Three  complete  steady  state  cycles  were  then  measured  for  input  voltages  ranging  from  600  V  to 
1600  V.  The  corresponding  field  inputs  to  the  model  were  determined  using  the  relation 

(4.1)  E  =  V/d 

where  d  =  0.0254  cm  is  the  thickness  of  the  wafer. 
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Table  3.  Coercive,  remanence  and  tip  characteristics  measured  from  200  mHz  data. 


PZT5A  (Circular) 
(1600  V) 

PZT5A  (Rectangular) 
(1600  V) 

PZT5H 

(2200  V) 

PZT4 

(1800  V) 

Ec  (V/m) 

1.2  x  10® 

1.3  x  10® 

0.87  x  10® 

1.4  x  10® 

Xc  ( C/(mV )) 

7.3  x  10-7 

7.2  x  10~7 

7.5  x  10"7 

1.2  x  10~® 

Pr  (C/m2) 

0.38 

0.38 

0.28 

0.31 

Xr  (< C/(mV )) 

4.4  x  10~8 

5.7  x  10"8 

8.2  x  10"8 

4.2  x  10-8 

Em  (V/m) 

6.2  x  10® 

4.2  x  10® 

5.8  x  10® 

4.7  x  10® 

Pm  (C/m2) 

0.46 

0.43 

0.38 

0.39 

Xln  (C/(mV)) 

2.8  x  10"8 

3.1  x  10“8 

1.5  x  10“8 

2.5  x  10“8 

Xm  (C/(mV)) 

5.0  x  10~9 

6.7  x  10~9 

2.1  x  10~9 

9.5  x  10~9 

The  parameters  a,a,c,k  and  Ps  were  estimated  using  both  techniques  discussed  in  Section  3.  The 
employment  of  the  data  characteristics,  summarized  in  Table  3,  in  Algorithm  2  provided  the  first  set  of 
parameter  values  listed  in  Table  2.  We  note  that  the  asymptotic  relations  employed  in  this  algorithm  are 
more  accurate  near  saturation  which  motivated  the  use  of  the  1600  V  input  data.  The  second  set  of  values 
were  obtained  through  a  least  squares  fit  to  the  1600  V  data.  A  comparison  of  parameter  values  obtained 
using  the  two  techniques  reveals  a  close  match  between  the  values  for  a,  a  and  k  with  some  discrepancy  in 
the  values  of  c  due  to  limitations  in  determining  Xm  and  Xm- 

The  modeled  polarization  was  then  computed  using  the  Langevin  anhysteretic  expression  and  the  pa¬ 
rameter  values  determined  through  the  least  squares  fit  to  1600  V  data  for  peak  input  voltages  ranging  from 
600  V  to  1600  V.  This  model  behavior  is  compared  with  the  measured  data  in  Figure  6.  We  first  note  that 
the  model  fit  for  the  1600  V  input  is  very  accurate  since  this  data  was  used  to  determine  the  parameters. 
Furthermore,  a  comparison  between  the  model  behavior  at  600  V,  800  V  and  1000  V  indicates  that  through¬ 
out  the  range  of  drive  levels,  the  model  very  accurately  predicts  the  measured  polarization.  We  reiterate 
that  this  capability  of  the  model  for  predicting  the  polarization  at  various  drive  levels  is  due  to  the  fact  that 
it  is  based  on  energy  principles. 

We  chose  to  employ  the  1600  V  data  for  the  least  squares  algorithm  solely  to  maintain  consistency  with 
the  asymptotic  Algorithm  2.  A  comparison  with  initial  results  presented  in  [20]  illustrates  that  for  this 
sample,  equally  accurate  model  fits  and  predictions  at  all  drive  levels  can  be  obtained  using  the  parameter 
values  a  —  3.7  x  106  Vm/C,  a  =  4.1  x  105  C/m2,  c  =  0.35,  k  =  1.8  x  106  C/m2  and  P8  =  0.49  C/m2  obtained 
through  a  least  squares  fit  to  the  600  V  data.  This  again  indicates  the  flexibility  of  the  model  due  to  its 
energy  formulation. 

Actuator  2  (Rectangular  Patch): 

The  second  PZT5A  actuator  which  we  considered  was  a  rectangular  patch  having  a  thickness  of  0.0381  cm 
(15  mils).  We  include  this  configuration  to  illustrate  the  variability  which  can  occur  between  batches  and 
the  manner  through  which  it  affects  the  model.  The  data  also  illustrates  certain  frequency- dependencies 
which  must  be  accommodated  in  broadband  applications. 
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To  illustrate  the  potential  variability  among  actuators,  the  data  collected  from  the  rectangular  patch 
with  a  peak  input  of  1800  V  is  compared  in  Figure  7a  with  1200  V  data  from  the  10  mil  thick  circular  patch. 
While  the  field  relation  (4.1)  indicates  that  the  polarization  should  theoretically  agree,  the  data  illustrates 
a  significant  difference  in  the  coercive  field  due  to  differences  between  the  materials.  This  necessitates  the 
refinement  of  parameters  to  attain  accurate  model  fits  throughout  the  range  of  operation. 

The  effect  of  even  slight  frequency  shifts  is  illustrated  in  Figure  7b  where  quasistatic  data  collected  at 
200  mHz  and  a  peak  input  voltage  of  2200  V  is  compared  with  corresponding  1  Hz  data.  (Note  that  this 
is  the  only  figure  where  1  Hz  data  is  included.)  In  its  current  state,  the  model  parameters  would  require 
updating  to  accommodate  the  observed  difference  in  saturation  behavior.  While  this  can  be  done  if  the 
frequency  is  fixed,  updating  in  this  manner  relies  on  the  mathematical  rather  than  physical  properties  of  the 
model.  The  extension  of  the  model  to  physically  incorporate  frequency  effects  is  under  investigation. 

The  modeled  polarization  curves  obtained  with  the  updated  parameters  summarized  in  Table  2  are 
compared  with  the  data  in  Figure  8.  It  is  observed  that  as  with  the  circular  patch,  the  model  accurately 
characterizes  the  hysteresis  throughout  the  drive  range  of  the  material. 


Example  2  (PZT5H): 

A  second  commonly  employed  soft  PZT  material  is  PZT5H.  In  this  example,  we  illustrate  the  per¬ 
formance  of  the  model  for  a  variety  of  parameter  choices  through  a  comparison  with  experimental  data 
generated  with  input  voltages  of  600  V,  1000  V,  1600  V  and  2200  V- 

The  initial  parameter  values  summarized  in  Table  2  were  obtained  using  Algorithm  2  and  the  least 
squares  algorithm  with  the  2200  V  data.  A  comparison  of  the  parameter  values  for  k  and  c  with  those  of 
PZT5A  indicates  smaller  values  for  the  PZT5H  sample.  The  respective  decreases  in  both  the  energy  required 
to  break  pinning  sites  and  the  amount  of  reversible  domain  wall  bending  reflects  the  fact  that  PZT5H  is 
softer  than  PZT5A  as  evidenced  by  the  lower  coercive  field.  We  also  note  that  for  this  soft  material,  the 
parameter  k  is  nearly  equal  to  the  coercivity  Ec  as  predicted  by  (3.15). 


Fig.  7.  Comparison  between  1200  V  data  from  the  0.0254  cm  (10  mil)  circular  patch  and  1800  V  data 
from  the  0.0381  cm  (15  mil)  rectangular  patch;  (b)  Comparison  between  1  Hz  and  200  mHz  data  from  the 
rectangular  patch. 
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Fig.  8.  Model  fit  to  600 ,  800,  1000,  1600  V  data,  collected  at  200  mHz  from  the  rectangular  PZT5A  patch, 
with  the  parameter  choices  c  =  0.15,  k  =  1.5  x  106  C/rr? ,  oc  =  3.7  x  106  Vm/G,  a  =  4.1  x  105  C/rr?  and 
Ps  =  .49  C/m2. 


The  model  fits  with  the  parameters  obtained  with  the  2200  V  data  are  compared  with  the  data  in 
Figure  9.  It  is  observed  that  the  model  with  parameters  estimated  using  the  least  squares  algorithm  provides 
an  excellent  characterization  throughout  the  hysteresis  cycle  whereas  the  model  with  parameters  specified 
by  Algorithm  2  accurately  quantifies  the  slope  at  the  remanence,  coercive  and  saturation  points  but  exhibits 
a  slight  discrepancy  in  polarization.  This  indicates  the  benefit  of  employing  a  least  squares  algorithm  to 
obtain  final  parameter  values. 

The  model  with  parameters  obtained  through  a  least  squares  fit  to  2200  V  data  is  then  used  to  predict 
the  polarization  at  lower  drive  levels  with  the  resulting  fits  plotted  in  Figure  10.  It  is  observed  that  with  these 
parameters,  the  model  provides  accurate  predictions  down  through  1000  V  but  degrades  somewhat  at  600  V. 
For  applications  which  require  the  full  range  of  actuator  operation,  a  second  strategy  is  to  consider  the  least 
squares  fits  to  data  from  a  variety  of  drive  levels  to  provide  parameters  which  optimize  the  model  performance 
throughout  the  input  range.  This  resulted  in  a  modification  of  the  value  of  k  from  k  =  1.0  x  106  C/m2  to 
k  =  0.9  x  106  C/m2  and  produced  the  model  fits  plotted  in  Figure  11.  With  this  choice,  the  characterization 
at  2200  V  is  slightly  less  accurate  but  the  model  prediction  at  600  V  is  improved.  We  reiterate  that  one 
set  of  parameters  is  still  employed  throughout  the  entire  input  range.  In  this  case,  however,  the  parameters 
have  been  optimized  for  the  full  operational  range  of  the  material. 

Finally,  we  note  that  the  parameter  values  predicted  by  Algorithm  2  using  the  2200  V  data  are  close 
to  the  values  obtained  through  a  least  squares  fit  due  to  the  accuracy  of  the  asymptotic  algorithm  relations 
near  saturation.  To  indicate  the  degradation  in  the  accuracy  of  the  relations  which  can  occur  at  lower 
drive  levels,  we  note  using  the  1000  V  data,  the  algorithms  produced  the  parameters  a  =  3.6  x  106  Vm/C, 
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a  =  4.3  x  105  C/m2,  it  =  1.2  x  106  C/m2  and  c  =  0.08  which  differ  significantly  from  the  values  obtained 
near  saturation  or  with  the  least  squares  fit.  This  indicates  an  advantage  of  the  least  squares  method,  which 
is  independent  of  drive  level,  and  the  necessity  of  applying  the  asymptotic  relations  near  saturation. 


Fig.  9.  Model  fits  to  the  PZT5H  data  with  parameters  obtained  through  the  least  squares  fit  and 
Algorithm  2;  (a)  Least  Squares:  a  —  4.2  x  106  Vm/C,  a  =  6.4  x  105  (7/m2,  k  =  1.0  x  106  C/m 2 ,  c  =  0.2, 
Ps  =  .425  C/m2,  (b)  Algorithm  2:  a  =  4.0  x  106  Vm/C,  a  =  5.8  x  105  C/m2,  k  =  1.1  x  106  C/m2, 
c  =  0.14, P8  =  .425  C/m2. 


Fig.  10.  Model  fit  to  the  PZT5H  data  with  the  parameter  choices  c  =  0.2,  k  =  1.0  x  106  C/m2, 
a  =  4.2  x  106  Vm/C,  a  =  6.4  x  105  C/m2  and  Ps  =  .425  C/m2. 
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Fig.  11.  Model  fit  to  the  PZT5H  data  with  the  parameter  choices  c  =  0.2,  k  —  0.9  x  106  (7/m2, 
a  —  4.2  x  106  Vm/C,  a  —  6.4  x  105  C/m2  and  Ps  =  .425  C/m2. 


Example  3  (PZT4): 

The  final  material  that  we  consider  is  the  harder  compound  PZT4  which  is  initially  poled.  Because 
more  energy  is  required  to  turn  dipoles,  the  material  is  not  depoled  when  cycled,  even  at  high  input  voltage 
levels.  The  remanent  bias  produces  the  asymmetry  observed  in  Figure  3  and  necessitates  the  inclusion  of 
the  bias  field  E0  and  polarization  values  P0  and  Pi  as  indicated  in  (2.6)  -  (2.8).  The  bias  can,  however,  be 
reduced  if  the  material  is  maintained  at  a  high  voltage  under  thermally  controlled  conditions.  The  subsequent 
trajectories  are  nearly  symmetric  and  exhibit  minimal  bias  in  field  or  polarization. 

In  this  example,  we  illustrate  the  capability  of  the  model  to  quantify  both  the  biased  (asymmetric) 
and  unbiased  (symmetric)  polarization  curves  for  PZT4.  We  consider  first  the  characterization  of  biased 
data  collected  at  a  peak  input  voltage  of  2200  V.  The  data  and  model  fit  obtained  with  the  parameter 
values  a  —  6.4  x  106  Vm/C,  a  —  8.0  x  105  C/m2,  k  =  1.5  x  106  C/m2,  c  =  0.5,  Ps  =  .44  C/m2  and 
E0  =  -4.0  x  105  V/m,  P0  =  .02  C/m2,  Pi  =  0  are  plotted  in  Figure  12.  It  is  observed  that  while  the 
model  is  not  able  to  completely  quantify  domain  switching  after  positive  remanence,  it  does  characterize  the 
primary  behavior  of  the  material  throughout  the  cycle  including  the  biased  field  and  polarization  and  the 
associated  asymmetry.  The  model  with  the  same  hysteresis  parameters  and  Eq  =  Pq  =  Pi  is  then  compared 
to  the  nearly  symmetric  data,  obtained  after  exposure  to  high  input  fields,  in  Figure  13.  There  it  accurately 
predicts  the  polarization  at  1000  V,  1200  V  and  1800  V  input  levels  but  under-predicts  the  hysteresis  present 
at  600  V.  We  note  that  these  fits  can  be  improved  if  the  parameters  are  refined. 

This  example  further  illustrates  the  philosophy  employed  in  this  model.  Through  the  incorporation 
of  the  associated  physics,  the  model  provides  a  characterization  dependent  only  on  inputs  to  the  actuator 
and  parameters  quantifying  the  state  of  the  material.  Once  the  hysteresis  parameters  a,o,c,  k  and  Ps  have 
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been  determined,  the  polarization  obtained  under  quasistatic  and  isothermal  operating  conditions  can  be 
characterized  using  the  measured  input  field  and  bias  values  Eq  and  P0.  This  provides  the  model  with 
significant  flexibility  in  a  variety  of  applications. 


Electric  Field  (MV/m) 


Fig.  12.  Model  fit  to  2200  V  PZT4  data  with  the  parameters  a  =  6.4  x  106  Vm/C ,  a  =  8.0  x  105  (7/m2, 
k  =  1.5  x  106  (7/m2,  c  =  0.5,  Ps  =  .44  <7/m2  and  E0  =  -4.0  x  105  V/m,  P0  =  .02  (7/m2,  Px  =  0. 


Fig.  13.  Model  fit  to  PZT4  data  with  the  parameters  a  =  6.4  x  106  Vm/C,  a  =  8.0  x  105  C/m2, 
k  =  1.5  x  106  (7/m2,  c  =  0.5,  Ps  =  .44  C/m 2  and  P0  =  Po  =  0. 
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5.  Concluding  Remarks.  This  paper  addresses  the  modeling  of  hysteresis  in  piezoelectric  materials 
through  the  application  and  extension  of  a  domain  wall  theory  for  ferroelectric  materials  [18,  19].  The  theory 
characterizes  the  inherent  hysteresis  in  the  relation  between  the  input  field  and  output  polarization  through 
the  quantification  of  energy  required  to  bend  and  translate  domain  walls  pinned  at  inclusions  in  the  material. 
This  provides  reversible  and  irreversible  polarization  components  whose  sum  represents  the  polarization  due 
to  an  applied  field.  Characterization  in  this  manner  provides  the  model  with  the  capability  for  specifying 
the  polarization  at  a  variety  of  input  field  levels  with  one  set  of  model  parameters.  The  flexibility  of  the 
model  is  further  augmented  by  the  small  number  (five)  of  required  parameters  and  the  physical  nature  of  the 
parameters.  For  example,  the  saturation  polarization  Ps  is  often  known  a  priori  or  can  be  directly  obtained 
from  the  data,  the  reversible  coefficient  c  can  be  estimated  from  the  ratio  of  the  slopes  of  the  polarization 
curve  at  field  reversal,  and  for  soft  materials,  the  pinning  coefficient  k  can  be  directly  estimated  from  the 
coercive  field  Ec- 

The  contributions  of  this  paper  are  threefold.  The  first  focussed  on  the  extension  of  the  model  to 
accommodate  hysteresis  loops  which  are  not  rotationally  symmetric.  This  was  motivated  by  the  biasing 
effects  observed  in  poled  hard  PZT  actuators  but  the  theory  is  sufficiently  general  to  encompass  a  variety 
of  applications.  The  second  contribution  was  the  validation  of  the  theory  and  illustration  of  its  predictive 
capabilities  for  three  commonly  employed  PZT  compounds.  The  theory  in  [18,  19]  was  illustrated  only 
for  PMN-PT-BT  employed  in  the  ferroelectric  regime  with  no  demonstration  of  the  predictive  capabilities, 
so  this  significantly  extends  the  applicability  of  the  theory.  The  third  contributions  of  the  paper  was  the 
development  of  asymptotic  relations  which  highlight  the  physical  nature  of  the  five  parameters  and  can  be 
used  to  obtain  initial  parameter  estimates.  While  the  approach  is  analogous  to  that  employed  in  [13]  for 
ferromagnetic  materials,  the  resulting  relations  and  algorithms  differ  in  certain  respects. 

The  model,  with  parameters  estimated  using  both  the  asymptotic  relations  and  a  least  squares  fit 
to  the  measured  data,  was  used  to  characterize  PZT5A,  PZT5H  and  PZT4  actuators.  In  each  case,  the 
parameters  were  determined  using  data  from  one  drive  level.  The  model,  with  these  parameters  fixed,  was 
then  used  to  predict  the  polarization  throughout  the  range  of  operation.  As  illustrated  in  the  examples,  the 
formulation  of  the  model  in  terms  of  energy  relations  based  on  the  input  field  provided  it  with  the  capability 
for  prediction  throughout  the  operational  range.  Furthermore,  once  the  biasing  field  Eo  and  polarization  P0 
were  determined,  the  model  could  accommodate  the  asymmetry  exhibited  by  the  poled  PZT4.  This  provides 
the  model  with  significant  flexibility  in  a  variety  of  applications. 

In  its  current  formulation,  the  theory  is  limited  to  quasistatic  and  thermally  controlled  operating  regimes. 
Moreover,  it  was  developed  under  the  assumption  that  the  effects  of  crystalline  anisotropies  are  minimal.  We 
note  that  in  certain  cases,  parameters  can  be  determined  which  provide  accurate  model  fits  at  a  variety  of  fixed 
frequencies.  The  use  of  the  model  in  this  manner  relies  on  its  mathematical  rather  than  physical  properties, 
however,  which  limits  its  flexibility  and  robustness  with  regard  to  changing  dynamics.  The  extensions  of  the 
physical  theory  to  accommodate  the  effects  of  frequency,  transient  temperatures  and  crystalline  anisotropies 
are  under  current  investigation. 

Finally,  the  ODE  nature  of  the  model  makes  it  amenable  to  inversion  through  the  consideration  of  a 
complementary  ODE  in  a  manner  analogous  to  that  described  in  [17].  This  facilitates  the  construction  of  an 
inverse  compensator  which  can  be  used  for  linear  control  design  [21,  22].  The  application  of  these  techniques 
for  linear  control  implementation  for  piezoceramic  actuators  which  exhibit  hysteresis  is  under  investigation. 
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